1 Setup

setwd("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/all-batches")
suppressPackageStartupMessages({
    library(data.table)
    library(DESeq2)
    library(gplots)
    library(here)
    library(hyperSpec)
    library(limma)
    library(LSD)
    library(magrittr)
    library(matrixStats)
    library(parallel)
    library(pander)
    library(plotly)
    library(RColorBrewer)
    library(scatterplot3d)
    library(tidyverse)
    library(tximport)
    library(VennDiagram)
    library(vsn)
})
suppressPackageStartupMessages({
    source("~/Git/UPSCb/UPSCb-common/src/R/featureSelection.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/gopher.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/plot.multidensity.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/volcanoPlot.R")
})
lfc <- 1
FDR <- 0.01
pal <- c(brewer.pal(8,"Dark2"),1)
pal2 <- brewer.pal(9,"Paired") #require package RColorBrewer
cols <- rainbow(17)
hpal <- colorRampPalette(c("blue","white","red"))(100)
mar <- par("mar")

2 Raw data

2.1 Loading of sample information

  • Read the sample information
samples <- read.csv("~/Git/UPSCb/projects/arabidopsis-nutrition-TOR/doc/samples3.csv")
  • Remove unnecessary samples
samples %<>% filter(!grepl("P11554_1",SciLifeID)) %>% 
    filter(! SciLifeID %in% c("P13406_101",
                        "P14066_128",
                        "P14066_133",
                        "P13406_102",
                        "P14066_131")) %>%
    mutate(Nutrition,Nutrition=relevel(Nutrition,"NPS")) %>% 
    mutate(AZD,AZD=relevel(AZD,"DMSO"))

samples <- samples[order(samples$Timepoint, samples$Nutrition, samples$AZD),]

samples %<>% mutate(Timepoint,Timepoint=factor(paste0("T",Timepoint)))

2.1.1 Load the data

  • Call the data
filenames <- list.files("../Salmon", 
                    recursive = TRUE, 
                    pattern = "quant.sf",
                    full.names = TRUE)
  • Name the data
names(filenames) <- sub("_S.*","",sapply(strsplit(filenames, "/"), .subset, 3))
  • Match data <=> sample list
filenames <- filenames[match(samples$SciLifeID,names(filenames))]
filenames <-filenames[!is.na(names(filenames))]
samples <- samples[match(names(filenames),samples$SciLifeID),]
  • Annotate the samples
samples$Conditions <- factor(paste(samples$Timepoint,samples$Nutrition,samples$AZD,sep="_"))
samples$Batch <- factor(substr(samples$SciLifeID,1,8))

3 Expression data

Read the expression at the transcript level

tx <- suppressMessages(tximport(files = filenames, 
                                type = "salmon", 
                                txOut = TRUE))

summarise to genes

tx2gene <- data.frame(TXID=rownames(tx$counts),
                      GENEID=sub("\\.[0-9]+","",rownames(tx$counts)))
gx <- summarizeToGene(tx,tx2gene=tx2gene)
## summarizing abundance
## summarizing counts
## summarizing length
kg <- round(gx$counts) 

Sanity check

stopifnot(all(colnames(kg) == samples$SciLifeID))

3.1 Raw data export

#dir.create(file.path("analysis_Tom","salmon"),showWarnings=FALSE,recursive=TRUE)
#write.table(kg,file="/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/nutrition-unormalised-gene-expression_data.csv")
#save(kg, samples, file = "/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/counts.rda")

3.2 Preliminary validations

3.2.1 Check for the genes that are never expressed

sel <- rowSums(kg) == 0 
sprintf("%s%% (%s) of %s genes are not expressed",
        round(sum(sel) * 100/ nrow(kg),digits=1),
        sum(sel),
        nrow(kg))
## [1] "16.9% (5452) of 32310 genes are not expressed"

4 Data normalisation

For visualization, the data is submitted to a variance stabilization transformation using DESeq2. The dispersion is estimated independently of the sample type

dds <- DESeqDataSetFromMatrix(
    countData = kg,
    colData = samples,
    design = ~ Conditions)
## converting counts to integer mode
dds <- estimateSizeFactors(dds)

4.1 Perform a Variance Stabilizing Transformation for plotting

vst <- varianceStabilizingTransformation(dds,blind=FALSE)
vsd <- assay(vst)
vsd <- vsd - min(vsd)

Write out

#write.csv(vst,"/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/analysis_Tom/library-size-normalized_variance-stabilized_data_nutrition.csv")

5 Multivariate analysis

5.1 PCA

Principal Component Analysis on the normalized data * Establishment of the PCA

conditions1 <- factor(paste(samples$AZD,samples$Nutrition,sep="_"))
pc <- prcomp(t(vsd))
percent <- round(summary(pc)$importance[2,]*100);percent
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 
##   46   21   14    6    3    2    1    1    1    0    0    0    0    0    0    0 
## PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29 PC30 PC31 PC32 
##    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
## PC33 PC34 PC35 PC36 PC37 PC38 PC39 PC40 PC41 PC42 PC43 PC44 PC45 PC46 PC47 PC48 
##    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
## PC49 PC50 PC51 PC52 PC53 PC54 PC55 PC56 
##    0    0    0    0    0    0    0    0
  • Graphical representation of PC1 x PC2
conds <- droplevels(conditions1)
plot(pc$x[,1],
     pc$x[,2],
     xlab=paste("Comp. 1 (",percent[1],"%)",sep=""),
     ylab=paste("Comp. 2 (",percent[2],"%)",sep=""),
     col=pal[as.integer(conds)],
     pch=c(15,16,17)[as.factor(samples$Timepoint)],
     main="All timepoints")

legend("bottomright",pch=19,
       col=pal[1:nlevels(conds)],
       legend=levels(conds))

legend("topleft",pch=c(15,16,17),
       col="black",
       legend=c("T0","T6","T24"))

  • Graphical representation of PC2 x PC3
plot(pc$x[,2],
     pc$x[,3],
     xlab=paste("Comp. 2 (",percent[2],"%)",sep=""),
     ylab=paste("Comp. 3 (",percent[3],"%)",sep=""),
     col=pal[as.integer(conds)],
     pch=c(15,16,17)[as.factor(samples$Timepoint)],
     main="All timepoints")

legend("topleft",pch=19,
       col=pal[1:nlevels(conds)],
       legend=levels(conds))

legend("bottomleft",pch=c(15,16,17),
       col="black",
       legend=c("T0","T6","T24"))

  • Graphical representation of PC1 x PC3
plot(pc$x[,1],
     pc$x[,3],
     xlab=paste("Comp. 1 (",percent[1],"%)",sep=""),
     ylab=paste("Comp. 3 (",percent[3],"%)",sep=""),
     col=pal[as.integer(conds)],
     pch=c(15,16,17)[as.factor(samples$Timepoint)],
     main="All timepoints")

legend("bottomleft",pch=19,
       col=pal[1:nlevels(conds)],
       legend=levels(conds))

legend("bottomright",pch=c(15,16,17),
       col="black",
       legend=c("T0","T6","T24"))

6 Differential expression based on the nutrition and treatment at T6

6.1 Filtration of samples based on timepoint

sel <- samples$Timepoint %in% c("T6")
suppressMessages(dds <- DESeqDataSetFromMatrix(
    countData = kg[,sel],
    colData = samples[sel,],
    design = ~ Nutrition * AZD))

6.2 Differential expression analysis

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

6.3 Variance Stabilising Transformation

6.3.1 Perform a Variance Stabilizing Transformation for plotting

vst <- varianceStabilizingTransformation(dds,blind=FALSE)
vsd <- assay(vst)
vsd <- vsd - min(vsd)

6.4 Contrasts to NPS_DMSO

The contrast by default is the first one (not Intercept)

resultsNames(dds)
## [1] "Intercept"            "Nutrition_NP_vs_NPS"  "Nutrition_NS_vs_NPS" 
## [4] "Nutrition_PKS_vs_NPS" "AZD_AZD_vs_DMSO"      "NutritionNP.AZDAZD"  
## [7] "NutritionNS.AZDAZD"   "NutritionPKS.AZDAZD"

6.4.1 Nutrition effect of carbon starvation

6.4.1.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "Nutrition_NP_vs_NPS")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
SucEffect6hrs <- rownames(res[cutoffs,])
SucLow6hrs <- rownames(res[cutoff2,])
SucHigh6hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 510 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 233 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 277 genes that are repressed

6.4.1.2 MA plot

DESeq2::plotMA(res)

6.4.1.3 Volcano plot

volcanoPlot(res, lfc = 1)

6.4.1.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.2 Nutrition effect of phosphorus starvation

6.4.2.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "Nutrition_NS_vs_NPS")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
PiEffect6hrs <- rownames(res[cutoffs,])
PiLow6hrs <- rownames(res[cutoff2,])
PiHigh6hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 3284 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 1246 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 2038 genes that are repressed

6.4.2.2 MA plot

DESeq2::plotMA(res)

6.4.2.3 Volcano plot

volcanoPlot(res, lfc=1)

6.4.2.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.3 Nutrition effect of nitrogen starvation

6.4.3.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "Nutrition_PKS_vs_NPS")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
NitEffect6hrs <- rownames(res[cutoffs,])
NitLow6hrs <- rownames(res[cutoff2,])
NitHigh6hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 290 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 162 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 128 genes that are repressed

6.4.3.2 MA plot

DESeq2::plotMA(res)

6.4.3.3 Volcano plot

volcanoPlot(res, lfc=1)

6.4.3.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

6.4.4 Nutrition effect of AZD-8055

6.4.4.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "AZD_AZD_vs_DMSO")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
AzdEffect6hrs <- rownames(res[cutoffs,])
AzdLow6hrs <- rownames(res[cutoff2,])
AzdHigh6hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 2564 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 871 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 1693 genes that are repressed

6.4.4.2 MA plot

DESeq2::plotMA(res)

6.4.4.3 Volcano plot

volcanoPlot(res, lfc=1)

6.4.4.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

7 Differential expression based on the nutrition and treatment at T24

7.1 Filtration of samples based on timepoint

sel <- samples$Timepoint %in% c("T24")
suppressMessages(dds <- DESeqDataSetFromMatrix(
    countData = kg[,sel],
    colData = samples[sel,],
    design = ~ Nutrition * AZD))

7.2 Differential expression analysis

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

7.3 Variance Stabilising Transformation

7.3.1 Perform a Variance Stabilizing Transformation for plotting

vst <- varianceStabilizingTransformation(dds,blind=FALSE)
vsd <- assay(vst)
vsd <- vsd - min(vsd)

7.4 Contrasts to NPS_DMSO

The contrast by default is the first one (not Intercept)

resultsNames(dds)
## [1] "Intercept"            "Nutrition_NP_vs_NPS"  "Nutrition_NS_vs_NPS" 
## [4] "Nutrition_PKS_vs_NPS" "AZD_AZD_vs_DMSO"      "NutritionNP.AZDAZD"  
## [7] "NutritionNS.AZDAZD"   "NutritionPKS.AZDAZD"

7.4.1 Nutrition effect of carbon starvation

7.4.1.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "Nutrition_NP_vs_NPS")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
SucEffect24hrs <- rownames(res[cutoffs,])
SucLow24hrs <- rownames(res[cutoff2,])
SucHigh24hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 718 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 503 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 215 genes that are repressed

7.4.1.2 MA plot

DESeq2::plotMA(res)

7.4.1.3 Volcano plot

volcanoPlot(res, lfc=1)

7.4.1.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

7.4.2 Nutrition effect of phosphorus starvation

7.4.2.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "Nutrition_NS_vs_NPS")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
PiEffect24hrs <- rownames(res[cutoffs,])
PiLow24hrs <- rownames(res[cutoff2,])
PiHigh24hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 4124 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 2167 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 1957 genes that are repressed

7.4.2.2 MA plot

DESeq2::plotMA(res)

7.4.2.3 Volcano plot

volcanoPlot(res, lfc = 1)

7.4.2.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

7.4.3 Nutrition effect of nitrogen starvation

7.4.3.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "Nutrition_PKS_vs_NPS")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
NitEffect24hrs <- rownames(res[cutoffs,])
NitLow24hrs <- rownames(res[cutoff2,])
NitHigh24hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 123 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 53 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 70 genes that are repressed

7.4.3.2 MA plot

DESeq2::plotMA(res)

7.4.3.3 Volcano plot

volcanoPlot(res, lfc=1)

7.4.3.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

7.4.4 Nutrition effect of AZD-8055

7.4.4.1 Extraction of the results from the DESeq analysis

res <- results(dds,name = "AZD_AZD_vs_DMSO")

cutoffs <- abs(res$log2FoldChange) >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff1 <- res$log2FoldChange >= lfc & ! is.na(res$padj) & res$padj <= FDR
cutoff2 <- res$log2FoldChange <= -lfc & ! is.na(res$padj) & res$padj <= FDR
AzdEffect24hrs <- rownames(res[cutoffs,])
AzdLow24hrs <- rownames(res[cutoff2,])
AzdHigh24hrs <- rownames(res[cutoff1,])

message(sprintf("There are %s genes that are differentially expressed",sum(cutoffs)))
## There are 5747 genes that are differentially expressed
message(sprintf("There are %s genes that are induced",sum(cutoff1)))
## There are 2876 genes that are induced
message(sprintf("There are %s genes that are repressed",sum(cutoff2)))
## There are 2871 genes that are repressed

7.4.4.2 MA plot

DESeq2::plotMA(res)

7.4.4.3 Volcano plot

volcanoPlot(res, lfc=1)

7.4.4.4 Heatmap

hm <- heatmap.2(t(scale(t(vsd[cutoffs,]))),col=hpal,
                Colv=FALSE,dendrogram = "row",trace = "none",labRow = FALSE,
                distfun = pearson.dist, labCol = paste(samples$Nutrition,samples$AZD,sep="-")[sel],
                hclustfun = function(X){hclust(X,method="ward.D2")},margins = c(6,5))

8 Comparisons of GOI lists (Venn diagrams)

8.1 GOI at T6

  • All GOI
grid.draw(venn.diagram(list(AzdEffect6hrs, PiEffect6hrs, SucEffect6hrs, NitEffect6hrs),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Phos. Starv.","Suc. Starv.", "Nit. Starv.")))

  • Downregulated genes
grid.draw(venn.diagram(list(AzdLow6hrs, PiLow6hrs, SucLow6hrs, NitLow6hrs),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Phos. Starv.","Suc. Starv.", "Nit. Starv.")))

  • Upregulated genes
grid.draw(venn.diagram(list(AzdHigh6hrs, PiHigh6hrs, SucHigh6hrs, NitHigh6hrs),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Phos. Starv.","Suc. Starv.", "Nit. Starv.")))

8.2 GOI at T24

  • All GOI
grid.draw(venn.diagram(list(AzdEffect24hrs, PiEffect24hrs, SucEffect24hrs, NitEffect24hrs),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Phos. Starv.","Suc. Starv.", "Nit. Starv.")))

  • Downregulated genes
grid.draw(venn.diagram(list(AzdLow24hrs, PiLow24hrs, SucLow24hrs, NitLow24hrs),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Phos. Starv.","Suc. Starv.", "Nit. Starv.")))

  • Upregulated genes
grid.draw(venn.diagram(list(AzdHigh24hrs, PiHigh24hrs, SucHigh24hrs, NitHigh24hrs),
                       filename=NULL,
                       col=pal[1:4],
                       category.names=c("AZD","Phos. Starv.","Suc. Starv.", "Nit. Starv.")))

8.3 Between timepoints

8.3.1 In response to AZD treatment

  • Upregulated
grid.draw(venn.diagram(list(AzdHigh6hrs, AzdHigh24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

  • DOwnregulated
grid.draw(venn.diagram(list(AzdLow6hrs, AzdLow24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

8.3.2 In response to Phosphate starvation

  • Upregulated
grid.draw(venn.diagram(list(PiHigh6hrs, PiHigh24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

  • DOwnregulated
grid.draw(venn.diagram(list(PiLow6hrs, PiLow24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

8.3.3 In response to Sucrose starvation

  • Upregulated
grid.draw(venn.diagram(list(SucHigh6hrs, SucHigh24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

  • DOwnregulated
grid.draw(venn.diagram(list(SucLow6hrs, SucLow24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

8.3.4 In response to Nitrogen starvation

  • Upregulated
grid.draw(venn.diagram(list(NitHigh6hrs, NitHigh24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

  • DOwnregulated
grid.draw(venn.diagram(list(NitLow6hrs, NitLow24hrs),
                       filename=NULL,
                       col=pal[1:2],
                       category.names=c("6hrs","24hrs")))

9 Export list of DEGs

write.table(AzdHigh24hrs, file = "AzdHigh24hrs.txt", sep = "\t", row.names = FALSE)
write.table(AzdHigh6hrs, file = "AzdHigh6hrs.txt", sep = "\t", row.names = FALSE)
write.table(PiHigh24hrs, file = "PiHigh24hrs.txt", sep = "\t", row.names = FALSE)
write.table(PiHigh6hrs, file = "PiHigh6hrs.txt", sep = "\t", row.names = FALSE)
write.table(SucHigh24hrs, file = "SucHigh24hrs.txt", sep = "\t", row.names = FALSE)
write.table(SucHigh6hrs, file = "SucHigh6hrs.txt", sep = "\t", row.names = FALSE)

write.table(AzdLow24hrs, file = "AzdLow24hrs.txt", sep = "\t", row.names = FALSE)
write.table(AzdLow6hrs, file = "AzdLow6hrs.txt", sep = "\t", row.names = FALSE)
write.table(PiLow24hrs, file = "PiLow24hrs.txt", sep = "\t", row.names = FALSE)
write.table(PiLow6hrs, file = "PiLow6hrs.txt", sep = "\t", row.names = FALSE)
write.table(SucLow24hrs, file = "SucLow24hrs.txt", sep = "\t", row.names = FALSE)
write.table(SucLow6hrs, file = "SucLow6hrs.txt", sep = "\t", row.names = FALSE)

10 Session Info

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] vsn_3.54.0                  VennDiagram_1.6.20         
##  [3] futile.logger_1.4.3         tximport_1.14.0            
##  [5] forcats_0.4.0               stringr_1.4.0              
##  [7] dplyr_0.8.3                 purrr_0.3.3                
##  [9] readr_1.3.1                 tidyr_1.0.0                
## [11] tibble_2.1.3                tidyverse_1.3.0            
## [13] scatterplot3d_0.3-41        RColorBrewer_1.1-2         
## [15] plotly_4.9.1                pander_0.6.3               
## [17] magrittr_1.5                LSD_4.0-0                  
## [19] limma_3.42.0                hyperSpec_0.99-20180627    
## [21] ggplot2_3.2.1               lattice_0.20-38            
## [23] here_0.1                    gplots_3.0.1.1             
## [25] DESeq2_1.26.0               SummarizedExperiment_1.16.0
## [27] DelayedArray_0.12.0         BiocParallel_1.20.0        
## [29] matrixStats_0.55.0          Biobase_2.46.0             
## [31] GenomicRanges_1.38.0        GenomeInfoDb_1.22.0        
## [33] IRanges_2.20.1              S4Vectors_0.24.0           
## [35] BiocGenerics_0.32.0         data.table_1.12.6          
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       rprojroot_1.3-2        htmlTable_1.13.2      
##  [4] XVector_0.26.0         fs_1.3.1               base64enc_0.1-3       
##  [7] rstudioapi_0.10        affyio_1.56.0          bit64_0.9-7           
## [10] AnnotationDbi_1.48.0   lubridate_1.7.4        xml2_1.2.2            
## [13] splines_3.6.1          geneplotter_1.64.0     knitr_1.26            
## [16] zeallot_0.1.0          Formula_1.2-3          jsonlite_1.6          
## [19] broom_0.5.2            annotate_1.64.0        cluster_2.1.0         
## [22] dbplyr_1.4.2           BiocManager_1.30.10    compiler_3.6.1        
## [25] httr_1.4.1             backports_1.1.5        assertthat_0.2.1      
## [28] Matrix_1.2-18          lazyeval_0.2.2         cli_1.1.0             
## [31] formatR_1.7            acepack_1.4.1          htmltools_0.4.0       
## [34] tools_3.6.1            affy_1.64.0            gtable_0.3.0          
## [37] glue_1.3.1             GenomeInfoDbData_1.2.2 Rcpp_1.0.3            
## [40] cellranger_1.1.0       vctrs_0.2.0            preprocessCore_1.48.0 
## [43] gdata_2.18.0           nlme_3.1-142           xfun_0.11             
## [46] rvest_0.3.5            testthat_2.3.0         lifecycle_0.1.0       
## [49] gtools_3.8.1           XML_3.98-1.20          zlibbioc_1.32.0       
## [52] scales_1.1.0           hms_0.5.2              lambda.r_1.2.4        
## [55] yaml_2.2.0             memoise_1.1.0          gridExtra_2.3         
## [58] rpart_4.1-15           latticeExtra_0.6-28    stringi_1.4.3         
## [61] RSQLite_2.1.2          highr_0.8              genefilter_1.68.0     
## [64] checkmate_1.9.4        caTools_1.17.1.2       rlang_0.4.2           
## [67] pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14         
## [70] htmlwidgets_1.5.1      bit_1.1-14             tidyselect_0.2.5      
## [73] R6_2.4.1               generics_0.0.2         Hmisc_4.3-0           
## [76] DBI_1.0.0              pillar_1.4.2           haven_2.2.0           
## [79] foreign_0.8-72         withr_2.1.2            survival_3.1-7        
## [82] RCurl_1.95-4.12        nnet_7.3-12            modelr_0.1.5          
## [85] crayon_1.3.4           futile.options_1.0.1   KernSmooth_2.23-16    
## [88] rmarkdown_1.18         locfit_1.5-9.1         readxl_1.3.1          
## [91] blob_1.2.0             reprex_0.3.0           digest_0.6.23         
## [94] xtable_1.8-4           munsell_0.5.0          viridisLite_0.3.0